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ABSTRACT 

A scheme is developed which enables one to trace backwards in time the cosmic density 
and velocity fields, and to determine accurately the current-epoch velocity field from 
the current-epoch density field, or vice versa. The scheme implements the idea of 
Giavalisco et al. (1993) that the principle of least action should be used to formulate 
gravitational instability as a two-point boundary-value problem. We argue that the 
Eulerian formulation of the problem is to be preferred to the Lagrangian one, on 
grounds of computational simplicity, of ease of interfacing with observational data, 
and of internal consistency at early times. The scheme is successfully tested on an 
exact solution in one dimension, and on currently Gaussian fields in one and two 
dimensions. The application of the scheme to real observational data appears to be 
eminently feasible, though computationally costly. 
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1 INTRODUCTION 

A key problem in cosmology is to understand and predict the 
evolution of the cosmic over-density and peculiar-velocity 
fields. Two propositions concerning these fields are widely 
accepted: 

I. The evolution of the cosmic density and velocity fields 
on scales significantly larger than that of an individual 
galaxy has been dominated by gravity. 

II. The cosmic density and velocity fields started from 
fluctuations of negligible amplitude at early times. 

Proposition I enables us, in principle, to predict the cur- 
rently observed fields that would have emerged from any 
given initial conditions. Conversely, it also enables us to 
determine the values of the fields at early times that are 
compatible with given fields at the current epoch. It follows 
that acceptance of both Propositions I and II is tantamount 
to placing strong constraints on the values of the fields at the 
current epoch. It turns out that these constraints leave only 
one of the fields free: once the density field has been chosen, 
the velocity field follows from the constraints, and vice versa. 
As was first pointed out by Peebles (1989), Proposition II 
in effect makes cosmology into a two-point boundary-value 
problem: observations provide us with a boundary condi- 
tion at the current epoch, while Proposition II imposes a 
boundary condition at t = from well-motivated physical 
assumptions. 

Given that the current velocity and density fields are in- 
dependently measurable, the capability to predict one field 
from the other provides a valuable empirical test of cos- 
mogony. Moreover, if the two fields prove to be related in 
the way predicted by theory, it will be of interest to inves- 
tigate the nature of the initial conditions. In particular, 



one would like to know whether the initial density field was 
Gaussian, as the simplest theory of the quantum origin of 
the cosmic density field in an inflationary universe predicts 
(e.g. Mukhanov et al. 1992; Peebles 1993). 

Although these things are possible in principle, it is not 
at all clear how one should set about doing them in practice. 
One set of major complications arises from the fragmentary 
and error-prone nature of observational data. This paper 
is not concerned with such considerations, but concentrates 
instead on the dynamical problem of relating current-epoch 
data to vanishingly small initial density fluctuations. 

Traditionally this problem has been approached through 
linear perturbation theory. Following Lifshitz (1946), one 
derives a linear o.d.e. for the evolution of the amplitude 
of each Fourier component 8h of the fractional over-density 
field 8(t,x). The general solution of this equation is a su- 
perposition of growing and decaying modes. Proposition II 
demands that the coefficient of the decaying mode be set 
equal to zero, with the result that the most general field 
compatible with Propositions I and II may at any time be 
represented as a sum of growing modes, with one arbitrary 
coefficient per mode. These coefficients precisely represent 
the independent degrees of freedom left to the fields by 
Propositions I and II. In fact, in linear theory, the veloc- 
ity field is simply a well-determined multiple of the gradient 
of the Newtonian gravitational potential, so that the general 
velocity field may be expressed in terms of the arbitrary co- 
efficients that determine the density field. 

Unfortunately, at the current epoch the fields are sig- 
nificantly non-linear on scales smaller than 8/t -1 Mpc and it 
is essential to have a theoretical framework that is valid in 
the non-linear regime. Zel'dovich (1970) made an important 
first step towards such a framework with his approximation 
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x(t,q)=x(Q,q)-D(t)V$\ <0t9) , (1) 

where x(t) is the position at time t of the particle q, D(t) is 
the linear growth factor that emerges from Lifshitz' o.d.e., 
and $ is the Newtonian potential. [It is conventional to la- 
bel particles by their initial (Lagrangian) coordinates: one 
sets q = a:(0,<jf).] Numerical experiments (e.g. Efstathiou et 
al. 1985, Coles et al. 1993) have demonstrated that, in spite 
of its simplicity, the Zel'dovich approximation is remarkably 
accurate. Except for one-dimensional systems, it is only an 
approximation, however, and at late times it fails to mimic 
the formation of gravitational potential wells around col- 
lapsing structures because in it particles move in straight 
lines. Hence it underestimates the density contrast in re- 
gions where clustering is significant. Moreover, for a given 
velocity field it leads to two, mutually incompatible, esti- 
mates of the density field, one obtained from the continuity 
equation and the other from Euler's equation of conservation 
of momentum. 

Several papers have recently proposed refinements to 
the Zel'dovich approximation that aim to extend the util- 
ity of the latter further into the non-linear regime (Nusser & 
Dekel 1992, Gramann 1993a,b, Buchert & Ehlers 1993). Per- 
haps because Zel'dovich obtained his approximation through 
inspired intuition rather than a systematic perturbation pro- 
cedure, no unique refinement has emerged, although it has 
been shown that the Zel'dovich approximation belongs to a 
class of solutions in the lowest-order parameterization of per- 
turbed particle orbits (Buchert 1992; Croudace et al. 1994). 

Moutarde et al. (1991) enhanced the coordinate map 
(1) by writing 

x(t,q) = q - DV<f>\ q + D 2 S (2) (q), (2) 

and this has been used to good effect by Bouchet et al. (1992) 
and Gramann (1993b), amongst others. Unfortunately, the 
algebra involved in these computations is laborious even at 
second order, and the prospect of calculating quantities to 
higher orders is not alluring. (For a detailed derivation of 
second and third order Lagrangian perturbation theory, see 
Buchert & Ehlers 1993.) Furthermore, these various pertur- 
bative approaches differ in their formalism and to a consid- 
erable extent in their predictions, with the result that appeal 
has been made to ra-body simulations to resolve differences 
between them (e.g. Mancinelli et al. 1993). 

More importantly, in these models, a solution of the 
governing equations is constructed as a Taylor series in the 
growth factor D(t). The range of convergence of such a 
Taylor series is surely quite small because 8 is known to 
become singular when the first caustics form; that is, for 
D = Do, a thoroughly finite number, roughly equal to the 
inverse of the largest eigenvalue of the velocity deformation 
tensor. It would be surprising if the Taylor series for 8 did 
not already break down when D is quite a lot smaller than 
Do- Hence, one is discouraged from pursuing the approaches 
of Bouchet et al. (1992) and Gramann (1993b) to higher 
orders not only by the considerable labour involved, but 
also by the fear that the higher-order terms, once obtained, 
would not prove meaningful at the current epoch. 

The principle of least action is the natural device for 
solving a two-point boundary value problem such as that 



involved in seeking cosmic fields compatible with given ob- 
servations and Propositions I and II. Not only does the least- 
action principle provide a natural mechanism for constrain- 
ing the Universe to be smooth at early times, but it also 
ensures that expansions of the dynamical variables such as 
8 in powers of D are constructed not as Taylor series but 
as least-squares approximations to the true solutions (see 
below). Since any function can be accurately approximated 
by a polynomial of sufficient degree, although very few func- 
tions have Taylor series with large radii of convergence, ap- 
proximations to the dynamical variables obtained from the 
least-action principle are very much more likely to be useful 
than those obtained by a Taylor-series approach. 

The principle of least action is most easily formulated 
for a system with a finite number of degrees of freedom, such 
as a collection of point masses or rigid bodies. Indeed, by 
treating each galaxy as a point mass, Peebles (1989, 1990) 
used the least-action principle to reconstruct the motion of 
the galaxies of the Local Group, while Dunn & Laflamme 
(1993) modelled neighbouring galaxies as extended, rigid ob- 
jects to study the acquisition of angular momentum through 
tidal interactions, in reasonable agreement with the present 
observational estimates. 

But at early times galaxies do not exist, and it cannot 
be appropriate to model the matter content of the Universe 
as a homogeneous distribution of massive bodies, whether 
pointlike or extended; until galaxies form, the only realistic 
formulation is in terms of a continuous field. If one adopts 
a Lagrangian approach, this field is the difference 

*(t,q) = x(t,q)-q (3) 

between the current and original positions of mass elements. 
In the Eulerian approach, the fields are the over-density 
and peculiar-velocity fields discussed above. Giavalisco et 
al. (1993) have formulated the problem for both types of 
field theory in the context of the principle of least action. 

In principle, the Lagrangian approach is to be preferred 
because remains finite at all times, whereas 8 grows with- 
out limit at orbit-crossing time. Nevertheless, from the prac- 
tical point of view the Eulerian approach enjoys several im- 
portant advantages. First, its vector field v is by Kelvin's 
theorem expressible as the gradient of a potential a. Conse- 
quently, the Eulerian theory involves only two unrestricted 
scalar fields. Second, in the Eulerian theory the gravita- 
tional potential $ is related to 8 by Poisson's simple equa- 
tion, whereas the relationship between $ and is complex 
in the extreme. Third, the current values of the Eulerian 
fields are more readily deduced from sparse and error-prone 
observations than are the Lagrangian displacement vectors 
of galaxies, with the result that the Eulerian theory is 
easier to apply to observational data than the Lagrangian 
theory would be. Finally, in Eulerian theory one may read- 
ily study the tidal forces on any given comoving region by 
calculating the shear and tidal tensors, whereas this under- 
taking is difficult in the Lagrangian framework. 

For these reasons, we have chosen to implement the Eu- 
lerian scheme that was outlined by Giavalisco et al. (1993). 
Section 2 sets out the basic formalism. Section 3 applies 
this formalism to an exactly soluble one-dimensional model 
and examines the results for some particular realizations of 
Gaussian random fields. Section 4 presents results for two- 
dimensional fields. Section 5 explains the relation of the 
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present scheme to similar schemes for following the evolu- 
tion of the cosmic fields. Section 6 sums up and describes 
how the scheme might be applied to real data. 

2 EULERIAN RECONSTRUCTION 

Giavalisco et al. (1993) show that the standard equations 
governing the evolution of cosmological perturbations in d 
dimensions may be derived from the action 



S 



dt / d d xC, 



(4) 



where x are comoving coordinates and the Lagrangian den- 
sity C is given by 

C = l( 1 + S)v 2 -<fS + a\s+-V-[(l + S)v]}- J^ $|2 2 .(5) 

Here 8 and $ are the fractional over-density and associated 
Newtonian potential, v is the peculiar velocity field, a(t) is 
the expansion factor of the universe and Pb(t) is the unper- 
turbed density. Varying the action with respect to v and 
the potential $ we find 

v(t, x) = — Va(t, x), 
K ' a K >' (6) 

V 2 $(i,:c) = 4:TrG Pb a 2 6(t, x). 

The first of these equations tells us that a is the velocity 
potential, while the second is simply Poisson's equation. 
Varying the action (4) with respect to a and 8, and then 
eliminating v with the aid of the first of equations (6) yields 







+ -[(l + 8)Va], 



(7) 



= a + 



2a 2 



These are the continuity and Bernouilli equations. By dis- 
carding the quadratic terms in equations (7), we recover the 
equations of linear theory. 

When v is similarly eliminated between equations (5) 
and (6), C becomes 

|V$|2 (8) 



2a 2 " SirGpta 2 

This form demonstrates that a and 8 are mutually conjugate 
variables. 

The solutions to the field equations (7), are constrained 
by the mixed boundary conditions 



8(0,x) = and 8(to, x) = 8o(x), 



(9) 



where to is the present time and <5o (x) is a function that may, 
in principle, be determined observationally. We express the 
fields as Fourier integrals of the form 

co {' d 

<f>(t,x) = 47rG Pb a 2 J2Mt) I ^ 

n = 

oo f d d k 

v(t,x) = a^2g„(t) J J7^y v >>,n e 1 *"", 

71 = 



(10) 



a(t, x) = a 



CO ft j 

i^ 9n{t) J j^r 



Cth,n e 



where 

f n (t) = D(t)[D(t)-l] n , 

r m W = 0, 1, 2, . . 

g n (t) = D(t)[D(t)-l] , 
with D(ta) = 1. In view of (6), we have 
Vh,n = ika h , n , 



-8h,n/k 2 



(11) 



(12) 



so the coefficients 8h,n and au,n determine all four fields. 
Notice that 8h,o and cth,o are the Fourier transforms of the 
current-epoch fields rather than of the early-time fields as 
in conventional perturbative cosmology (e.g. Padmanabhan 
1993; Jain & Bertschinger 1994). Notice also that the 8h,n 
and cth,n do not depend on time, the time-dependence of 
the physical fields having been absorbed into the functions 
f n (i) and g n (i) defined by equations (11). 

Substituting the expansions (10) into equations (7), 
multiplying the equations through by g r and f r , respectively, 
and integrating over time, we obtain a nonlinear system of 
time-independent equations for the 8h,n and ah,n- 



~^2{fng r )8h,n - k 2 ^ ( 

71 = 71 = 

d d p 



71,771 = 



(Ug 



(13) 



(2*) 



k p8 k - 



■p,n p,m 



and, 



4irG > (pbfnfr)8k,n 



n = 



k 2 Y^((ffn + 2^-g n )f r )a k ,n 



n = 



E 



(14) 



d d p 
J2^Y 



p ■ (k -p) a h - P ,n a P , 



phere 



dt a (i)x(t) 



(15) 



The values of the angle brackets required in (13), and (14) 
are calculated in Appendix A for the case Q = 1, A = 0. The 
corresponding values for other cosmological models are eas- 
ily calculated. For concreteness, we shall henceforth assume 
that Q = 1, A = 0. 

Equations (13) and (14) provide a pair of equations for 
each wavevector k and non-negative integer r. Thus, for ev- 
ery k there is a infinite system of equations r = 0,1,2,.... 
We simplify this system by (i) limiting the index r to r < N , 
and (ii) assuming that the fields satisfy periodic boundary 
conditions in a box of side L. This assumption enables us 
to replace the Fourier integrals (10) by sums over a count- 
able set of wavenumbers k t . These sums we truncate after 
K wavevectors in each coordinate direction. With these re- 
strictions and approximations we have 



n = 



(16) 



a(t, x) 



a 2 (t 
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where the 8h,n and au,n satisfy (see Appendix B for inter- 
mediate results), 



N 

E 

n = 



<f>(n + r) 



2r - 3ra 
2(n + r) 



8k,n — k Cth,n 



5 N 

' oTd + m + ^ • P 6 h-p,n «p,m, 

P 



(17) 



2L d ^ 

n,m = 



and 

N 

^</>(ra + T) | Ofc,„ 



3r — 2ra 2 

K ttfc.n 



5^ 
6£ d 



3(ra + r) 

iV 

^ »?(« + m + r) • (fc - p) a h -p, n 



(18) 



Here 



0(n) = (-4) 



„ n\ (ra + 2)! 
(2ra + 5)! ' 
0(« + l) 



(19) 



n + 1 

For concreteness we restrict ourselves to the case in 
which the current-epoch density field is specified - the case 
of specified velocity potential is entirely analogous and re- 
quires only rearrangement of the coefficients. In this case 
it is expedient to arrange equations (17) and (18) such that 
the non-linear terms and terms involving 8h,o appear on the 
right. Then for each wavevector h we have the (2N + 1)- 
dimensional matrix equation 

A \ f fl{k) 
A T D^j \A(k) 



C (1 \k) 



(20) 



where the vectors to be calculated are, 

6fc,2 



fl{k) 



( a >>,o \ 

«fc,2 



A{k) 



\ »k,N / 

and the elements of the submatrices are 
A rn = <t>(r + n), 



(21) 



<f>(r + n), 



(22) 



(1) = 3r-2» 
u ™ ~ 3(« + r) 

2(n + r) 

The quadratic subcolumns on the right-hand side of (20) are 
given by, 

5k 2 N 

= ——^ ^ v( n + m + r) p ■ (k - p) a h - P , n « P ,m, 



<2) 



n,m = 
P 

N 

"2i7 E 

n,m = 
V 



(23) 



5 



r](n + m + r) k ■ p 8h- 



•p,n <-*p,m, 



while 

*r = <f>(r) 

independent of k. 



(24) 



The couple the matrix equations of the set (20) for 
different wavevectors k. So one would ideally solve (20) as 
a set of K d equations in K d (2N + l)-dimensional variables, 
using the Newton-Raphson algorithm to handle the non- 
linearity inherent in the C^'K For a realistic number K of 
wavevectors, this approach is unfeasible, however, and we 
have solved the system as follows. 

Given the Fourier coefficients 8h,o of a current-epoch 
density field, we estimate cth,o from standard linear theory: 



D, 



V 2 a 



k (<i 



(25) 



Then we set 8h,n = <Xh,n = for < n < N and evaluate 
the non-linear terms C^'K The r.h.s. of (20) may now be 
completely evaluated and fl(k), A(k) solved for by multi- 
plying (20) through by the inverse of the matrix on its l.h.s. 
The values of 8h,n and au,n for n ^ recovered in this way 
are used to re-evaluate the and the procedure iterated. 
One finds that the iterated values of the fields converge to 
a solution of the equations provided the initially specified 
fields are not large, and therefore the non-linear terms are 
of only secondary importance. 

If the initial fields are large, brute-force iteration of the 
type just described does not yield a solution, and we proceed 
as follows. We use brute-force iteration to solve the equa- 
tions obtained by replacing the by m~ 1 C^'\ where m is 
a sufficiently large integer. Then we use the solution we have 
just obtained as the starting point for an iterative solution 
of the equations obtained by replacing by 2m~ 1 C^'\ 

For m sufficiently large, these iterations also converge, en- 
abling the procedure to be repeated with 3m -1 on the 
r.h.s., and so on until equations (20) have been solved. 



3 APPLICATION TO ONE-DIMENSIONAL 
SYSTEMS 

3.1 An analytically soluble model 

Equations (20) simplify greatly in the case d = 1. More- 
over, in d = 1 we have an exact solution to the perturbative 
cosmological problem against which to test our apparatus. 
Namely the plane-wave solution (e.g. Padmanabhan 1993) 



8{t,x) 



(t/r) 2 ' 3 cos( K g) 



(26) 



1 - (t/r) 2 / 3 cos( K g) ' 

where the Eulerian coordinate x is related to the comoving 
Lagrangian coordinate q by 

' t \ 2 / 3 sin(Kg) 



x(t, q) 



(27) 



As t — > 0, the density field (26) tends to a sinusoidal wave. 
We investigate the ability of our machinery to recover this 
from the non-sinusoidal field described by (26) for t ~ r. 
At to the density contrast So given by (26) satisfies 



(toM 



2/3 



< 8 < 



(toM 



2/3 



(28) 



l + (Vr) 2 / 3 " " l-(Vr) 2 / 3 ' 

Hence, for t = r/2 3/2 ~ 0.35r, we have that So lies between 
<*>max = 1 and <5 m i n = —1/3. Evidently, nonlinear effects are 
important from at least this time onwards. 

The full curves in Fig. 1 show the analytical density 
field at t = (thin curve) and t = 0.46r (thick curve), while 
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Figure 1. Reconstruction of a one-dimensional density field. The 
current-epoch field (thick curve) is given by (26) with to = 0.46r 
and k = 27r/256. The order of the truncated series in (16) is 
N = 10. The thin curve shows the exact initial density from (26), 
while the dashed curve shows the reconstructed initial density. 
Notice that the evolution predicted by linear theory has been 
factored out. 

the dashed curve shows the initial density field obtained 
by applying the least-action principle to the field at time 
to = jr. Fig. 2 shows the corresponding velocity fields. 
Note that the quantities plotted are 8/D and v/aD, so any 
evolution apparent in these figures is due to non-linearity; 
the growth in amplitudes predicted by linear theory has been 
factored out. 

An accurate description of the sharp peaks and flat- 
bottomed troughs of 8(to,x) requires significant power in 
Sk,o for A large. When the reconstructed initial field is 
evaluated, this high-A power must be cancelled by an equal 
amount of high-A power in Sk,n for n > 0, to leave a sim- 
ple sinusoid. Failure to solve the system (20) accurately, 
either because one has not iterated to convergence or one 
has adopted too small a value of N, leads to incomplete 
subtraction of the high-A components of 8(to, x) and small- 
scale wiggles in the reconstructed fields. The latter are more 
noticeable in 8(0, x) than in v(0, x) since the high-A Fourier 
coefficients of v are smaller than those of 8 by a factor ~ A -1 . 

3.2 One-dimensional Gaussian fields 

Gaussian random fields are of special interest since they 
are the simplest random fields, and they arise naturally in 
the simplest models of the generation of cosmic fluctuations 
(e.g. Peebles 1993). They are characterized by a power spec- 
trum since the phases of their Fourier modes are uncorre- 
cted with one another. An important geometrical property 
of a Gaussian field 8(x) is that it is statistically indistin- 
guishable from —8(x). In particular, regions of given over- 
density have, on the average, the same physical shapes as 
regions of equivalent under-density. As an initially Gaussian 



Figure 2. The velocity fields associated with the density fields of 
Fig. 1. The full curves show the initial (thin) and current (thick) 
fields from (26), while the dashed curves show the corresponding 
fields obtained from the least-action principle. The current ve- 
locity field predicted by the least-action principle is plotted as a 
dotted curve, but is barely visible because it overlays the analytic 
curve. The velocity scale is such that the initial velocity ampli- 
tude is normalized to unity. The evolution predicted by linear 
theory has again been factored out. 

density field evolves under the action of gravity, it ceases to 
be a Gaussian field since over-dense and under-dense regions 
evolve differently. In A-space this is reflected in the develop- 
ment of correlations between the phases of the field's Fourier 
components. 

Soda et al. (1992) and Gramann (1992) have studied the 
growth of these correlations in numerical simulations. The 
interpretation of the results of these calculations is difficult, 
however, given that phases are determined only modulo 2ir. 

In this section we investigate the backwards evolution 
of fields that are Gaussian at the current epoch. In partic- 
ular, we ask whether a current-epoch Gaussian field can be 
reached by the action of gravity on a physically realizable 
initial field. We set 



Ah 



-2 2jrir (c 



(29) 



where < ru < 1 is a random number and the amplitude A 
is related to the variance in the density, a 2 , by 

42 



L 2 



A 2 L 2 
(2^ 



(30) 



C(4), 



where £ is Riemann's zeta function 

CO 

C(m) = ]T ra - m ; C(4)« 1.082. 



(31) 



The full curve in Fig. 3 shows a Gaussian field with 
Fourier coefficients given by (29), while the dashed curve 
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shows the corresponding initial field from equations (20). 
The current-epoch field selected shows a void at the centre, 
and a peak near x = (which is identified with x = 256). In 
the initial field, the void is narrower and the peak wider than 
at the current epoch. Relative to linear theory, the void is 
initially deeper, and the peak lower, than at present. Thus 
we are here simply seeing the usual tendency of gravity to 
widen voids and make them flatter-bottomed, whilst raising 
and sharpening peaks. These trends are apparent in all the 
realizations of random fields that we have studied. 

Fig. 4 shows the velocity fields associated with the den- 
sity fields of Fig. 3. The velocity changes sign near the 
extrema of 8, as we expect, given that velocity vectors point 
away from voids. Similarly, the largest absolute velocities 
are attained near points at which 8 passes through zero, 
which shift outward with the expansion of the void. 



4 TWO-DIMENSIONAL GAUSSIAN FIELDS 

The evolution of cosmic fields is much more interesting in 
two dimensions than in one, since two-dimensional fields dis- 
play phenomena such as shear and tidal fields that have 
no one-dimensional analogues (Ellis 1971; Matarrese, Pan- 
tano & Saez, et al. 1993; Bertschinger 1993). As Melott & 
Shandarin (1989) point out, two-dimensional fields have the 
advantage over three-dimensional ones of being readily dis- 
played and thus providing a useful laboratory in which to 
study the evolution of cosmic fields. 

We have again traced backwards currently Gaussian 
fields. We take 



Ak 



-3 27riru 

e h 



(32) 



where the phases ru are random numbers < < 1. A is 
related to the variance of the density field by 



E 



i 



A 2 L 

(2tt) 6 ^ U 2 + 

i,j=0 
i = j?0 



(33) 



: 4.66-" 



, A 2 L 2 
W 



The full contours in Fig. 5 are contours of constant den- 
sity at the current epoch of a field that is now Gaussian with 
dispersion <j = 1.1. A peak at center left is balanced by a 
broad, flat depression that is deepest near (20, 5). The con- 
tour 8 = is the heavy line. The reconstructed early-time 
field (dashed for 8 > 0, dotted for 8 < 0) shows the peak 
to have formed by the amalgamation of two density max- 
ima, one centred on (10,17) and a smaller one centred on 
(25, 16). At early times the trough contains two sharp min- 
ima at the base of a conical hole. Over time these minima 
merge and the hole becomes shallower and flatter-bottomed. 
Thus again we see the tendency of gravity to make troughs 
broader and shallower, and peaks taller and sharper, as well 
as its tendency to amalgamate similar features. 

Fig. 6 shows three velocity fields associated with the 
density fields of Fig. 5. The left-hand panel shows the pre- 
dicted current velocities. As expected, these diverge from 
the void and converge on the peak. The centre panel shows 
the reconstructed initial velocity field. This is similar to 
the current field, the main differences being (i) a shift in the 
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Figure 3. Backwards evolution of a current-epoch Gaussian field 
with Fourier coefficients given by (29) for a = 0.26. Full curve: 
current-epoch density field. Dashed curve: reconstructed initial 
field. The order of the truncated series (16) is N = 10. 
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Figure 4. Velocity fields corresponding to Fig. 3. Full curve: 
current-epoch velocity amplitude predicted by equations (20). 
Dashed curve: reconstructed initial velocity amplitude. 

main centre of convergence consequent on a shift in the loca- 
tion of the dominant peak, and (ii) a more complex structure 
consequent on the existence of a subsidiary peak and trough. 
The right panel in Fig. 6 shows the velocity field predicted 
by the Zel'dovich approximation - each arrow in the cen- 
tre panel has been displaced the appropriate amount, and 
where this places more than one arrow in a cell, the arrows 
in that cell have been averaged. Clearly, the velocity field 
shown in the right panel is a mediocre approximation to the 
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Figure 5. Contours of constant present (full) and initial (broken) 
density for a field which is now Gaussian with Fourier coefficients 
given by equation (32). The field's dispersion is (T = 1.1 and the 
lattice of Fourier modes has K = 32 points on a side. Contours 
are spaced by ^ and contours of negative initial over-density are 
dotted, while those of zero and positive over-density are dashed. 
The thick full curve marks the contour on which 8 = at the 
current epoch. 

actual current-epoch field shown in the left panel. In par- 
ticular, distortion by the subsidiary peak is still apparent, 
and the main point of convergence does not coincide very 
accurately with the actual density peak. In fact, the middle 
panel offers at least as good an approximation to the left 
panel, which suggests that the "frozen-flow" approximation 
(Matarrese et al. 1992) provides as accurate an estimate of 
evolved velocity fields as does the more widely employed 
Zel'dovich approximation. 

5 RELATIONSHIP TO OTHER WORK 

In this section we clarify the relation of the present approach 
to similar work in the literature. Several entirely different 



expansion schemes for following the evolution of the cosmic 
fields have been explored. Since these schemes often bear 
considerable superficial resemblance to one another, it may 
be felt helpful to summarize their relationship to this work. 

The most fundamental distinction is between Eulerian 
and Lagrangian schemes. Buchert (1992) and Buchert & 
Ehlers (1993) have explored schemes that are fully La- 
grangian in the sense that auxilliary fields are expressed as 
functions f(q) of the Lagrangian coordinate q. Moreover, 
this work differs from that of Bouchet et al. (1992), Gra- 
mann (1993b) and others in that q does not coincide with 
the Eulerian coordinate x at t = but at t = to > 0. Con- 
sequently, one can write 

x(t,q) =q + p(t,q), (34) 

where p — > as t — > to. The smallness of p is exploited 
to solve the non-linear equations of motion iteratively, with 
non-linear terms evaluated using p from the last iteration. 
The convergence of the series of solutions thus generated 
is by no means guaranteed. In fact, our own experience 
with the algebraically similar system (20) suggests that it 
will converge only for \t — to\ small. The solutions obtained 
diverge as t — > because they include decaying modes. Con- 
sequently, they are incompatible with Proposition II of the 
Introduction and cannot be used to investigate initial con- 
ditions. 

A large number of papers employ a similar iterative ap- 
proach within Eulerian theory (Juszkiewicz (1981); Goroff 
et al. 1986; Bernardeau 1992; Makino, Sasaki & Suto, 1992; 
Jain & Bertschinger 1994); the Fourier transforms of the 
Eulerian equations of motion and continuity are written so 
that linear terms appear on the left-hand sides and non- 
linear terms appear on the right. The fields 8h etc are then 
expanded as power series in D(t) and the resulting system of 
non-linear equations is solved iteratively with the right-hand 
sides evaluated in advance using the values of the fields from 
the previous iteration. The first approximation, being ob- 
tained by neglecting the right-hand sides, is simply standard 
linear theory with an arbitrarily chosen amplitude. As in the 
Lagrangian schemes described above, there is no guarantee 
that the series obtained converge, and little understanding 
of the optimum number of terms to employ. The key dif- 
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Figure 6. Velocity fields corresponding to the density field of Fig. 5. Left: the reconstructed current-epoch velocity field. Centre: the 
reconstructed initial velocity field. Right: the velocity field predicted by applying the Zel'dovich approximation to the reconstructed 
initial fields. 
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ferences between this approach and the present work are 
(i) that we obtain an exact solution of our truncated equa- 
tions rather than generating a probably asymptotic power 
series, and (ii) we express the early-time fields as functions 
of a current-epoch field, rather than expressing the current- 
epoch fields as functions of the early-time field 8(0, x). The 
latter dependence is convenient if one wishes to determine 
ensemble averages of the fields rather than the evolution 
of a particular field, and one is willing to assume that the 
early-time fields are Gaussian. Our approach is the one to 
choose if one wishes to work with actual, observationally 
determined, fields. 

As explained in the Introduction, our work differs from 
that of Bouchet et al. (1992), Gramann (1993a, b) and others 
in that we seek least-squares approximations to the evolu- 
tion equations, rather than Taylor-series developments of 
these solutions. 



6 CONCLUSIONS 

The growth of inhomogeneities from a near-homogeneous 
early Universe constitutes a two-point boundary- value prob- 
lem and as such, it is most naturally approached through 
the least-action principle. The Eulerian formulation of the 
problem is simpler than the Lagrangian formulation in that: 
(a) it can be formulated in terms of two unrestricted scalar 
fields; (b) in it the gravitational potential is simply related 
to one of these fundamental fields; (c) its predictions for the 
current epoch may be easily and directly compared with ob- 
servations. The principal weakness of the Eulerian approach 
is that the density field becomes singular as soon as caustics 
cross, with the result that the radius of convergence of any 
Taylor series for the density contrast is probably small. For- 
tunately, when the variational constraint equations are em- 
ployed, the power series involved are least-squares fitted to 
the solutions rather than established as Taylor series, with 
the result that their validity is not as narrowly restricted 
as when one approaches the problem from the differential 
equations of motion. 

In Section 2 we have formulated the application of the 
least-action principle to cosmic inhomogeneities in terms of 
coefficients 8h,n and au,n in expansions of the Fourier am- 
plitudes 8h(t) and ah(t) in powers of (D — 1) (n = 0, 1, . . .), 
where D is the linear growth factor and a(t, x) is the veloc- 
ity potential. These coefficients solve a non-linear system of 
equations (20). In these equations one may consider either 
the amplitudes 8h,o or the amplitudes cth,o to play the role 
of boundary conditions, which specify the structure of the 
cosmic fields at the present time. The solutions to the equa- 
tions then yield the structure of the conjugate field (a or 8, 
respectively) at the current epoch as well as the structure of 
both fields at all earlier epochs. 

In Section 3 we have solved the system (20) for the one- 
dimensional case, for which an exact solution is known. We 
obtain solutions by increasing the magnitude of the non- 
linear terms in the equations to their true values in stages, 
and iterating a trial solution to a fully converged solution 
at each stage. This technique enables us to solve equations 
(20) for any meaningful values of the current-epoch fields. 

The calculated backwards evolution of the exact so- 
lution agrees well with the form of the solution at earlier 
times provided sufficient (N > 3) powers of (D — 1) are 



employed. A field whose Fourier representation involves sig- 
nificant power at large wavenumbers k evolves backwards 
into a pure sinusoid in consequence of the high-A compo- 
nents of 8(to,x) being more and more exactly cancelled by 
the high-A coefficients of (D — l) n for large n. Errors in 
the completeness of this cancellation determine the exacti- 
tude with which early-time fields can be determined. The 
early-time and current-epoch velocity fields can be more ac- 
curately determined because they involve less high-A power. 

We have used our machinery to trace backwards in one 
and two dimensions fields that are currently Gaussian. We 
find that at early times their peaks are lower and flatter 
topped, and their troughs deeper and more pointed. In two 
dimensions this exercise reveals the tendency of structures 
to form by the merging of sub-structures. Hence, on trac- 
ing a peak backwards in time in Section 4, we found it to 
have formed by the merging of two smaller peaks that were 
initially located on either side of it, while a broad shallow 
trough had emerged from the amalgamation of two rather 
vee-shaped holes. 

Since our reconstruction scheme yields both the current- 
and the early-epoch velocity fields, it enables us to test the 
validity of the Zel'dovich approximation. In Section 4 we 
concluded that the Zel'dovich approximation does not pro- 
vide a better approximation than the "frozen-flow" approx- 
imation of Matarrese et al. (1992). 

Section 5 clarified the relation of this work to similar 
expansion schemes that have been employed to evolve the 
cosmic fields forwards in time. The essential differences are: 
(i) while we approximate solutions of the evolution equations 
by polynomials in (D — 1), other schemes generate power- 
series approximations that are either Taylor series or formal 
power series that are probably at best asymptotic; (ii) we ex- 
press the early-time fields as functions of the current-epoch 
fields, while in many other schemes the functional depen- 
dence is reversed. 

Before the reconstruction scheme developed here can 
be applied to real data, two things have to be done: (i) 
the scheme has to be implemented in d = 3, (ii) it has to 
be interfaced with an effective means of determining either 
8h,o or ah,o- Implementation in d = 3 is straightforward; 
the only drawbacks to working in three dimensions are the 
difficulty of displaying the results and the greatly increased 
computational cost - both the memory and the CPU time 
required are proportional to K d . In fact, we have already 
tested a three-dimensional implementation of the present 
scheme. 

The observational determination of either 8h,o or cth,o 
is harder. The aj^o are the easiest quantities to estimate, 
since every distance to a galaxy of known redshift yields an 
estimate of f g ■ Va, where f g is the direction to that galaxy. 
An attractive scheme is to choose the ah o which minimize 



E 



k 



(35) 



where x g and v g are the position and radial peculiar velocity 
of the c/th galaxy in some sample. Since this is a standard 
linear least-squares problem, its solution would be straight- 
forward if one's galaxy sample contained more galaxies than 
a reasonable grid in fc-space does points (> 30, 000). In prac- 
tice, one is unlikely to be able to observationally determine 
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so many coefficients cth,o in the near future, and it will be 
necessary to set aj^o = for k > fc max , say. Nonetheless, for 
n > one's calculations should include coefficients au,n for 
k > fc max , since as one traces the fields backwards in time, 
power will appear at these high frequencies even if none is 
present in a(x) at the current epoch. We hope to report 
on an investigation along these lines in the not too distant 
future. 
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APPENDIX A: EVALUATION OF ANGLE 
BRACKETS 

We evaluate for the special case Q = 1 the integrals over 
time required in equations (13) and (14). By the definition 
(15) of (x), they are 



(fn.gr 



dt a(tf fn(t) g r (t), 



{g n gr) = I dt a(t) g n {t) g r {t), 
Jo 

(fng r g m ) = dt a(tf f n {t) g r {t) g m {t), (Al) 
Jo 

\2 



(Pbfnfr 



dt a(ty Pb (t) f n (t) f r (t), 



((gn + 2-g n )f r )= J dta(tf [g n (t) + 2-g n (t)]f r (t), 
where 



\ 2/3 

«tt=l£) and p b (t) = 



(A2) 



We use the result D = a = (t/to) 2 ^ 3 to effect the replace- 
ment 



dta -+ |*o / dD D° 



5/2 . (A3) 
o Jo 

Next use the definitions (11) to evaluate /„, g n and their 
time derivatives in terms of D . This yields 

2 



{fn9r) = 3^ [ 1{ l' n + T) + nI & n + r " 1} 
2 

(g n gr) = n + r )> 



3tn 



{fng m g r ) = T7--f(f , m + n + r), 



3t 



(A4) 



{Pbfnfr) 



-J(f,n + r), 



4wGt ^ 2 ' 

((g n + 2% n )f r ) = ±- [/(!,« + r) + |n7(|, 
tt to 

where for n = 0, 1, 2, ... and a > 



n + r — 1) 



I(a, n) = I du u a (u — l) n 



(A5) 



(The proof of the last equality is by induction on n from 
n = by repeated integration by parts.) It follows that 



/(f,n-l) = -A/(!,n). 
Since 

n+l 

TJg I a r i { , 2 2n (n + 2)! 
1H»+ 2 J " 4 - 2 ( 2 „ + 5)!' 



(A6) 



(A7) 
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we have finally 



7(|, n) = 4!0(ra) , J(§ , n - 1) = — 0(ra), (A8) 



phere 



0(n) = (-4)"" !(w + 2 >, ! . (A9) 
vy ' v 7 (2ra + 5)! v 7 



APPENDIX B: CONSTRAINT EQUATIONS 

On substituting from Appendix A for the values of the an- 
gle brackets, truncating the sums over ra, and replacing the 
Fourier integrals by Fourier sums as described in the text, 
equations (13) and (14) become 



N 

E 

n = 



<f>(n + r) 



1 - ^ 

1 2 



ra 



k Cth,n 



ra + r i 

E 5 <j>(m + ra + r + 1) 
2 ra + ra + r + 1 

-^j ^ • p 8 h - p ,n «p,m, 



iV 

^ </>(« + r) 

n = 



5ra 



3(ra + r) 



E 



5A; 2 0(m + ra + r + 1) 
6 m + ra + r + 1 



(Bl) 



X ~P) -P «fc-p,m a P ,». 



The sum over p in the first of equations (Bl) can be 
rewritten as a convolution since 

^^kp 8 h - p a p = ^ [(k-p)-p 8 h - p a p + 8 h - p a p p 2 ] .(B2) 



When we confine the Fourier sums to a sublattice of K 
points, these convolutions can be evaluated efficiently, and 
without further approximation, by application of the dis- 
crete Fourier transform theorem. 
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